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Abstract. We investigate properties of dense suspensions and sediments of small spherical silt 
particles by means of a combined Molecular Dynamics (MD) and Stochastic Rotation Dynamics 
(SRD) simulation. We include van der Waals and effective electrostatic interactions between the 
colloidal particles, as well as Brownian motion and hydrodynamic interactions which are calculated 
in the SRD-part. We present the simulation technique and first results. We have measured velocity 
distributions, diffusion coefficients, sedimentation velocity, spatial correlation functions and we have 
explored the phase diagram depending on the parameters of the potentials and on the volume 
fraction. 
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I. INTRODUCTION 

We simulate clay like colloids, for which in 
many cases the attractive Van-der- Waals forces 
are relevant. They are often called "peloids" 
(Greek: clay- like). The colloidal particles have 
diameters in the range of some nm up to some 
/im. In general, colloid science is a large 
field, where many books have been published 
Ull H II HI! The term "peloid" originally 
comes from soil mechanics, but particles of this 
size are also important in many engineering pro- 
cesses. Our model system of ALjOs-particles of 
diameter 0.5 /xm suspended in water is an of- 
ten used ceramics and plays an important role 
in technical processes. In soil mechanics0 and 
ceramics science Q, questions on the shear vis- 
cosity and compressibility as well as on porosity 
of the microscopic structure which is formed by 
the particles, arise0,0|. In both areas, usu- 
ally high volume fractions ($ > 20%) are of 
interest. The mechanical properties of these 
suspensions are difficult to understand. Apart 
from the attractive forces, electrostatic repul- 
sion strongly determines the properties of the 
suspension. Depending on the surface poten- 
tial, which can be adjusted by the pH-value of 
the solvent, one can either observe formation 
of clusters or the particles are stabilized in sus- 
pension and do sediment only very slowly. Hy- 
drodynamic effects are also important for sed- 
imentation experiments. Since typical Peclet 
numbers are of order one in our system, Brow- 
nian motion cannot be neglected. 

In summary, there are many important fac- 
tors which have to be included into a model 
which describes peloids in a satisfying way. 
Such a model is needed to gain a deeper un- 
derstanding of the dynamics of dense colloidal 
suspensions. A lot of effort has been invested by 
applying different simulation methods, which 



have their inherent strengths but also some 
disadvantages. Simplified Brownian Dynamics 
(BD), such as in the work of H utter [TT| does 
not contain long-ranged hydrodynamic interac- 
tions among particles at all. The computational 
cost is low, since hydrodynamics is reduced to 
a simple Stokes force and thus large particle 
numbers can be handled. BD with full hydro- 
dynamic interactions utilizes a mobility matrix 
which is based on the Oseen- or Rotne-Prager- 
Yamakawa tensor approximations which are ex- 
act in the limit of zero Reynolds number and 
zero particle volume fraction 1 1 2lll3|. 
This technique faces the main problem that the 
computational effort scales with the cube of the 
particle number due to the inversion of matri- 
ces. 

The lattice Boltzmann (LB) method on the 
other hand is numerically efficient and in- 
trinsically contains hydrodynamic interactions. 
Ladd and Verberg give an overview over the 
LB method and describe how to include stress 
fluctuations 01- Adhikari et al. add noise to 
their model by introducing a noise term for ev- 
ery lattice velocity and node |15|. However, the 
discussion about the correct inclusion of ther- 
mal fluctuations is still ongoing [TEl Ht| . Pair- 
Drag simulations have been proposed by Silbert 
et al.|l7j. which include hydrodynamic interac- 
tions in an approximative way. They have fo- 
cused on suspensions with high densities (up 
to 50 %) of uncharged spherical colloidal parti- 
cles. Here we use Stochastic Rotation Dynam- 
ics (SRD) 0, , a recently developed method 
to simulate fluid flow, and combine this with 
a Molecular Dynamics (MD) simulation for the 
suspended particles. SRD is a particle-based 
method which does not show any numerical 
instabilities, contains thermal fluctuations in- 
trinsically and is simple to implement. Many 
important issues in fluctuation fluid dynamics 
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such as sedimentation [2(j , vesicles in flow |2j , 
polymers in flow [2*3, reacting fluids [23] have 
been addressed very recently using this method. 
In this paper, first we discuss the main points 
of the MD-part of our simulation code, sec- 
ond we present the SRD method in the context 
of our work, then we describe two alternative 
ways of coupling the two parts of the simula- 
tion and point out the advantages and disad- 
vantages of these two possibilities. After that, 
we analyze the time scales which are relevant for 
the peloids, we want to simulate. Based on the 
insights of this section we show in the following 
section how to determine the simulation param- 
eters. Then we describe how we have tested our 
simulation code and present first results in the 
following section. Finally in the last section we 
draw a conclusion and summarize shortly the 
model we have presented. 



II. MOLECULAR DYNAMICS 




center - center distance / particle radius 

FIG. 1: DLVO Potentials for AI2O3 spheres of 
R — 0.5 /im diameter suspended in water. These 
are typical potentials used for our simulations as de- 
scribed below. The primary minimum at d/R — 2.0 
is not reproduced correctly by the DLVO theory. It 
has to be modeled separately. In most of our cases 
the existence of the secondary minimum determines 
the properties of the simulated system. 



The colloidal particles in our simulation are 
represented by three dimensional spheres. In 
order to correctly model the statics and dy- 
namics when approaching stationary states, re- 
alistic potentials are needed. The interaction 
between the par ticles is described by DLVO 
theory [ToL ITll .l2"ij . If the colloidal particles are 
suspended in a solvent, typically water, ions 
move into solution, whereas their counter ions 
remain in the particle due to a different resolv- 
ability. Thus, the colloidal particle carries a 
charge. The ions in solution are attracted by 
the charge on the particles and form the electric 
double layer. It has been shown (see |23), that 
the resulting electrostatic interaction between 
two of these particles can be described by an 
exponentially screened Coulomb potential 
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where d denotes the particle diameter and r is 
the distance of the particle centers, z is the 
charge of the ions, e the elementary charge, 
T the temperature, \&n denotes the effective 
surface potential, and k is the inverse Dcbyc 
screening length. In addition the behavior is 
determined by the attractive van der Waals in- 
teraction which can analytically be integrated 
over the two spheres. This leads to the second 
part of the DLVO potential: 
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Ah is the Hamaker constant which involves the 
polarizability of the particles and of the solvent. 
The DLVO potentials are plotted in Fig.Q] for 
six typical examples with different depth of the 
secondary minimum. The primary minimum 
has to be modeled separately, as discussed be- 
low. 

To avoid that the particles penetrate each 
other, one needs a repulsive force depending on 
their overlap. We are using a Hertz force de- 
scribed by the potential 

Vkortz = K(d - rf/ 2 if r<d, (3) 

where K could be expressed by the elastic mod- 
ulus of AI2O3. This would determine the simu- 
lation time step, but to keep computational ef- 
fort relatively small, we determine the time step 
using the DLVO-potentials as described later 
on and then choose a value for K. Two aspects 
have to be considered: K has to be big enough 
so that the particles do not penetrate each other 
by more than approximately 10% and it may 
not be too big, so that numerical errors are kept 
small, which is the case when the collision time 
is resolved with about 20 time steps. Otherwise 
total energy and momentum are not conserved 
very well in the collision. 

Since DLVO theory contains the assumption of 
linear polarizability, it holds only for large dis- 
tances, i.e. the singularity when the two spheres 
touch, does not exist in reality. Nevertheless, 
there is an energy minimum about 30 fceT deep, 
so that particles which come that close would 
very rarely become free again. To obtain nu- 
merical stability of our simulation, we model 
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this minimum by a parabolic potential, some 
ksT deep (e.g. 6ksT). The depth of the mini- 
mum in our model is much less than in reality, 
but the probability for particles to be trapped 
in the minimum has to be kept low enough so 
that only few of them might escape during sim- 
ulation time. 

Long range hydrodynamic interaction is taken 
into account in a separate simulation for the 
fluid as described below. This can only repro- 
duce interactions correctly down to a certain 
level. On shorter distances, a lubrication force 
has to be introduced explicitly in the molecular 
dynamics simulation as described in 25]. The 
most dominant mode, the so-called squeezing 
mode, is an additional force 

F lub - -(v rel ,f)f 67rT??Ved , (4) 
r - ri - r 2 

with r rcd = ri7 " 2 (5) 
n + r 2 

between two spheres with radii r±, r 2 and the 
relative velocity v re i. r\ is the dynamic viscos- 
ity of the fluid. Fi u b diverges if particles touch 
each other. Therefore, we limit the force by in- 
troducing a minimum radius, where the force 
reaches its largest allowed value. The potential 
is shifted accordingly to smaller particle dis- 
tances, so that the maximum force is reached 
for particles touching each other. 
The Hertz force also contains a damping term 
in normal direction, 

FDamp = -(V r el, ?)f/fyr ~ ?1 - T 2 , (6) 

with a damping constant (3 and for the trans- 
verse direction a viscous friction proportional 
to the relative velocity of the particle surfaces 
is applied. 

For the integration of the translational mo- 
tion we utilize a velocity Verlet algorithm [2^ 
chap. 3.2.1 to update the velocity and position 
of particle i according to the equations 

Xi(i + 5t) = Xl (t) + St^(t) + St 2 ^-, (7) 

m 

2m 

For the rotation, a simple Euler algorithm is 
applied: 

u)i{t + 8t) = Ui(t)+StTi, (9) 
9i(t + 5t) = 9i(t) + F(9i,Wi,5t), (10) 

where u)i(t) is the angular velocity of particle 
i at time t, Tj is the torque exerted by non 
central forces on the particle i, 9i(t) is the ori- 
entation of particle i at time t, expressed by a 



quaternion, and F(0i,uji, St) gives the evolution 
of 9i of particle i rotating with the angular ve- 
locity u>i(t) at time t. 

The concept of quaternions [2(| is often used to 
calculate rotational motions in simulations, be- 
cause the Euler angles and rotation matrices 
can easily be derived from quaternions. Using 
Euler angles to describe the orientation would 
give rise to singularities for the two orientations 
with 9 = ±90°. The numerical problems re- 
lated to this fact and the relatively high com- 
putational effort of a matrix inversion can be 
avoided using quaternions. 

We have switched off dissipative forces and 
checked if the total energy and each component 
of the total momentum are conserved. We have 
verified this for the molecular dynamics simula- 
tion for the simulation of the fluid, and for the 
coupled simulation separately. 
We also checked that our implementation of the 
molecular dynamics code is correct by simulat- 
ing eight large particles with Hertz-repulsion 
and Coulomb friction in a closed box at a vol- 
ume fraction of <& ~ 20%. We checked that the 
collisions are realistic, i.e. that the individual 
angular velocities for two particles interacting 
in a non-central collision before and after they 
have touched are consistent. 



III. STOCHASTIC ROTATION 
DYNAMICS (SRD): SIMULATION OF 
THE FLUID 

The Stochastic Rotation Dynamics method 
(SR D) introduced by Malevanets and Kapral 
[LsL [I3 is a promising tool for a coarse-grained 
description of a fluctuating solvent, in par- 
ticular for colloidal and polymer suspensions. 
The method is also known as "Real-coded Lat- 
tice Gas" [23 or as "multi-particle-collision dy- 
namics" (MPCD) 28]. It can be seen as a 
"hydrodynamic heat bath", whose details are 
not fully resolved but which provides the cor- 
rect hydrodynamic interaction among embed- 
ded particles [29j. SRD is especially well suited 
for flow problems with Peclet numbers of or- 
der one and Reynolds numbers on the parti- 
cle scale between 0.05 and 20 for ensembles 
of many particles 01 ■ The numerical effort 
scales linearly with the number of embedded 
colloidal particles unlike in Brownian Dynam- 
ics, and only one random number per node (for 
the choice of the rotation matrix) is needed in 
contrast to fluctuating lattice-Boltzmann. For 
Peclet-numbers of order one, about three to five 
SRD-particles are required per box (or node) 
whose positions and velocities can be seen as 
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the degrees of freedom in that node. In three 
dimensions this amounts to 18 to 25 variables 
per node which is similar to the 15 or 19 speed 
lattice-Boltzmann method. 
While the LB method might be slower than 
SRD in the regime of large thermal fluctu- 
ations it has the advantage that it can be 
used for almost arbitrarily high Pcclct-numbers 
just by reducing the amplitude of the noise. 
To reduce the noise in SRD, a huge num- 
ber of fluid-particles per node has to be used 
which makes the method inapplicable at Peclet- 
numbers higher than about 20. Fortunately it 
has been shown by Louis and Padding 20] , that 
basic properties of sedimentation such as the 
main settling speed are hardly affected by ther- 
mal noise. 

The method is based on so-called fluid par- 
ticles with continuous positions and velocities 
which follow a simple, artificial dynamics. 

The system is coarse-grained into cubic cells 
of a regular lattice with no restriction on the 
number of particles in a cell. The evolution of 
the system consists of two steps: streaming and 
collision. In the streaming step, the coordinate 
of each particle is incremented by its displace- 
ment during the time step. Collisions are mod- 
eled by a simultaneous stochastic rotation of 
the relative velocities of every particle in each 
cell. The dynamics is explicitly constructed to 
conserve mass, momentum, and energy, and the 
collision process is the simplest consistent with 
these conservation laws. It has been shown that 
there is an H— theorem for the dynamics and 
that this procedure yields the correct hydrody- 
namic equations for an ideal gas|l8|. 

Consider a set of N point-particles with (con- 
tinuous) coordinates and velocities Vj(i). 
In the streaming step, all particles are propa- 
gated simultaneously by a distance v,r, where 
r is the value of the discretized time step. For 
the collision step, particles are sorted into cells, 
and they interact only with members of their 
own cell. Typically, the simplest cell construc- 
tion consisting of a hyper cubic grid with mesh 
size a, is used. The collision step consists of an 
independent random rotation of the relative ve- 
locities v— u, of the particles in each cell, where 
the macroscopic velocity u(£, t) is the mean ve- 
locity of the particles in the cell with coordinate 
£. The local temperature T(£,i) is defined via 
the mean square deviation of the particle ve- 
locities from the mean velocity in a cell. All 
particles in a cell are subject to the same ro- 
tation, but the rotation angles of different cells 
are statistically independent. There is a great 
deal of freedom in how the rotation step is im- 
plemented 0,H3, since, by construction, the 



local momentum and kinetic energy are invari- 
ant. The dynamics is therefore summarized by 



Ti(t + T) = Vi(t)+TVi(t), (11) 

Vt(t + r) = u[£i(t + T)]+w[Zi(t + T)](12) 

•{Vi(t)-U[&(t+T)]}, 



where u>(£i) denotes a stochastic rotation ma- 
trix, and & is the coordinate of the cell occu- 
pied by particle i at the time of the collision. 
u(£) = v/j is the mean velocity of the par- 

fee? 

tides in cell £. u> is taken to be a rotation by an 
angle ±a, with probability 1/2. We are using 
rotations about the three coordinate axes with 
a = ±90°, because these are the most simple 
rotation matrices one can imagine in 3D, since 
they only contain entries taken out of {0, zfclj. 
This has been suggested by M. Straufi in |31j . 
In every time step for each cell one of these 
6 possibilities is chosen with equal probability 
1/6. However, any stochastic rotation matrix 
consistent with detailed balance can be used. 

In order to remove low temperature anoma- 
lies and to achieve exact Galilean-invariance, we 
use a modification of the original algorithm |30j : 
all particles are shifted by the same random vec- 
tor with components in the interval [—a/2, a/2] 
before the collision step. Particles are then 
shifted back by the same amount after the col- 
lision. The random vectors of consecutive it- 
erations are uncorrelated. Ihle and Kroll have 
discussed in Ref. 32, 33) why this simple proce- 
dure works and shown that it leads to transport 
coefficients independent of an imposed homoge- 
neous flow field. In Ref. [34( and [35| analytical 
calculations of the transport coefficient in this 
method are presented. 

Two different methods to couple the SRD 
and the MD simulation have been introduced 
in the literature. We have implemented them 
both and we are using them depending on what 
we plan to measure. The first one|27| is much 
more accurate in resolving the local velocity 
field around the colloidal particles. Lubrica- 
tion effects are reproduced well by this coupling 
method. The second one[3^ resolves the veloc- 
ity field only down to a length scale of the parti- 
cle diameter. On the other hand the method be- 
comes much faster because of the lower resolu- 
tion. In both coupling methods the long range 
hydrodynamic interactions are reproduced. 
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IV. COUPLING I: PLACING FLUID 
PARTICLES OUTSIDE OF COLLOIDAL 
PARTICLES 

In the combined MD and SRD simulation the 
fluid particles have to interact with the colloidal 
particles and transfer momentum from one to 
the other part of the simulation. One possibil- 
ity to do this is, as suggested by Inoue et al.[27j 
to check after each streaming step of a fluid par- 
ticle i, if its new position Xi(t + r) is within a 
colloidal particle and if yes, to modify its po- 
sition and velocity. In this coupling step total 
momentum has to be conserved, which means, 
that when modeling the "collision" between the 
fluid particle and the colloidal particles, one has 
to make sure that the change of momentum of 
the fluid particle is transfered to the suspended 
particle. The calculations described in the fol- 
lowing are done in a frame fixed on the colloid 
particle. 

One can think of several different methods 
to assign a new position to the fluid particle, 
which have been shown to work properly: 

1) place it on the shortest distance to the sur- 
face of the colloidal particle and move it with 
its new velocity half of a time step, 

2) calculate the point and the exact time when 
the fluid particle has entered into the colloidal 
particle and move it back to there. Then choose 
a new velocity and move the fluid particle with 
the new velocity for the remainder of the time 
step. 

Both methods turned out to work, where the 
second one is more accurate but more compu- 
tationally intensive as well. Just to place the 
fluid particle directly on the surface and move 
it again in the next time step turned out to 
produce an increase of the fluid particle den- 
sity around the colloidal particle. Anomalies in 
the fluid temperature could also be found when 
the fluid particles were placed directly on the 
colloid surface. 

To increase stability of the simulation the 
idea is not to conserve energy in every single col- 
lision, but to use a thermostat and choose the 
new velocities according to a given distribution. 
The new velocities should point from the colloid 
surface to the outer area. Since the interior of 
the colloidal particle usually does not contain 
any fluid particles and the velocity distribution 
next to a colloidal particle should be indepen- 
dent of neighboring particles, the velocity dis- 
tribution for the newly chosen fluid particle ve- 
locities has to be the same as if the space in- 
side the suspended particle was filled with fluid 
particles. Assume these imaginary fluid parti- 
cles having the same density and temperature 



as in the remainder of the fluid bath. Then, 
one could evaluate the velocity distribution for 
the reflected fluid particles by taking the veloc- 
ity distribution of the imaginary fluid particles 
passing through the colloid surface. But it is 
a non-trivial task to analytically calculate this 
distribution for a spherical area. However, if 
the mean free path of the fluid particles is small 
compared to the diameter of the colloidal parti- 
cles, we can safely assume the colloid surface to 
be an infinitely extended plane separating the 
space into two regions j27|. Then one finds the 
following distribution: 

p{v n ) ~ u„exp(-/3^), (13) 
p{v t ) ~ exp(-/?t£) (14) 



where v n is the normal component and v t is 
the tangential component of the fluid particle 
velocity in the frame fixed to the surface of the 
large particle, m/ is the mass of a fluid parti- 
cle. In the following sections we describe how 
rrif has to be chosen. T is the temperature 
to which this thermostat is adjusted and the 
whole system will adopt this temperature af- 
ter a transient time. The tangential component 
can be obtained by computing \J x\ + x\ of two 
independent and Gaussian distributed random 
variables. 

Since the fluid particles of the SRD are arti- 
ficial particles within the context of this meso- 
scopic simulation method, their mean free path 
and their momentum are different from the cor- 
responding values for single solvent molecules. 
Because of this, there is a depletion force act- 
ing on colloidal particles which is much larger 
than in reality. Depletion forces are only rel- 
evant in systems with very big molecules, e.g. 
polymer solutions with added small particles or 
binary mixtures of particles with clearly sepa- 
rated diameters. There, each of the small par- 
ticles carries a considerable momentum - which 
is also the case in the SRD simulations. Never- 
theless, unrealistically high depletion forces can 
be suppressed by reflecting fluid particles many 
times: If after the collision step the fluid par- 
ticle is placed in another colloidal particle, the 
collision step is repeated for that colloidal par- 
ticle and so on, until the fluid particle reaches 
a position outside any colloidal particle or un- 
til a maximum number iV max of collisions has 
been calculated through. We have measured 
the depletion force and found out that a limit of 
Af max ~ 10 is a good compromise between com- 
putational speed and accuracy. The depletion 
force does not decay substantially stronger if 
the limit is increased, but the computational ef- 
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fort still grows with iV max (at most linearly), be- 
cause some fluid particles are trapped in a small 
gap between two colloidal particles and jump 
from one to the other. This in fact would still 
decrease the depletion force, but in the mean 
time the calculation for the remaining system 
is interrupted until finally eventually one single 
fluid particle is reflected the very last time. It 
is obvious that this scenario can easily be trun- 
cated. The remaining depletion force can be 
neglected at least in the cases where strong at- 
tractive van der Waals forces or strong repulsive 
electrostatic forces are present. 

V. COUPLING II: ROTATING 
VELOCITIES OF THE COLLOIDAL 

PARTICLES 

A second possibility to couple the SRD and the 
MD simulation is to sort the colloidal particles 
into the SRD boxes and include their velocity 
in the rotation step. This technique has been 
used to model protein chains suspended in a 
liquid f36l l37j ]. The mean velocity in each cell 
has then to be weighted with the mass of the 
particle (because the mass of colloid particles 
differs at least by one order of magnitude from 
the one of the fluid particles and their inertia 
dominates the flow field next to it). The cal- 
culation of u[£j(£ + t)] in Eq. Ijl2|l is modified 
to 

u (0 = Jj^2 v km k , (15) 

where we sum over all colloid and fluid parti- 
cles in the cell, is the mass of the particle 
with index k and M = ^ is the total mass 

contained in the cell. 

The coupling acts on the center of mass of 
the colloidal particles and affects only the fluid 
particles within the same cell. This means, to 
affect the same area of the flow field like in real- 
ity, one has to choose the cells to be of the same 
size as the colloidal particles. Obviously, the 
mesh size is drastically larger than in the first 
coupling method and the flow field cannot be 
resolved in detail. The fact that colloidal parti- 
cles push away the solvent as well as depletion 
and lubrication forces cannot be reproduced at 
any level. 

VI. TIME SCALE ANALYSIS 

Our system contains many different, let us say 
L, time scales, which differ by several orders of 



magnitude making brute force numerical simu- 
lations very time-consuming or even impossible. 
These time scales can be used to define L — 1 
dimensionless characteristic numbers, such as 
the Reynolds- or the Peclet number as the ra- 
tio of two time scales. If one can manage to 
adjust the simulation parameters such that all 
these characteristic numbers are the same as in 
the experiment, the simulations should be able 
to exactly reproduce the dynamical behavior of 
the real system. Of course, therefore one has to 
change quantities like the temperature or the 
viscosity of the fluid. 

Often, it is sufficient to reproduce only a few 
of all characteristic numbers exactly, i.e. only 
those which are believed to be significant for 
the behavior. For example, in sedimentation 
processes where the Reynolds number is much 
smaller than unity, it may be modified to an- 
other value, which still fulfills the condition of 
being much smaller than one. In both cases 
the Stokes limit is a valid approximation. As a 
general rule of thumb, dimensionless numbers 
of order one are important to be reproduced 
since they represent two competing dynamical 
effects. The reason to modify the other "in- 
significant" numbers is to reduce the ratio of the 
largest to the smallest time scale which deter- 
mines the numerical effort. In order to decide 
which are the dimensionless numbers that can 
be safely modified without changing the physics 
too much, a detailed analysis of the different 
time scales is needed. 

We start with the largest scales. After some 
time an isolated spherical particle sedimenting 
in a liquid reaches the so-called Stokes velocity, 

„_§*»(&_!). (16) 

9 v \p w J 

v is the kinematic viscosity, g denotes grav- 
ity! Pm is the mass density of the particle, p w 
the mass density of the solvent. This veloc- 
ity is obtained from the force balance between 
buoyancy and weight of the particle, Fa = 
47r(/0 m — p w )gR 3 /3 and the drag- force in a vis- 
cous liquid, Fjj — Qitvp w Rv. 
The drag-force Fr> also defines the mobility 
p = v/Fd = l/(6TTvp w R) of a spherical parti- 
cle. The time for a particle to move a distance 
of its diameter, 2R, is denoted by 

2R §v 
T s = — = 7 r- (17) 

By means of the Einstein-relation D — pk^T 
we obtain the diffusion constant D for the par- 
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tide, 



D = 



k B T 



(18) 



The mean square displacement of a diffusing 
particle in each dimension % is given by 



< x 2 {t) >=2Dt, 



(19) 



hence, the time the particle needs to diffuse a 
distance of 2R is of the order of 



td 



2R 2 



12nvp w R 3 



(20) 



which we call diffusion time. 
The ratio Tp, /t$ measures the importance of 
Brownian motion in the system and is called 
Peclet number, Pe = — . It turns out to be 
close to unity here. Inserting the definitions for 
td and ts, one notices that Pe depends on the 
fourth power of the radius R, 



Pe = 



v s R _ F G R _ ^gR\ Pn 



Pv 



D 



3k B T 



(21) 



Let us consider another time proportional to 
td : we assume a regular three-dimensional, cu- 
bic array of spheres which are separated by gaps 
of size R/2. Then, the volume concentration of 
this suspension is in the intermediate regime, 
4> = 0.268. The time one sphere diffuses the 
distance of a gap is given by tq = td/16. 
Now, let us discuss another important time 
called the particle relaxation time, which is re- 
lated to how long it takes the particle to react 
to an imposed force, i.e. this time measures the 
inertia effects. Consider Newton's equation for 
a particle of mass m subject to a force F and a 
friction coefficient £, = —t;v + F. Expand- 
ing the velocity v around the stationary state, 
v = vs + 5v, gives 



dSv £ . 

~»r = — dv ' 

at m 



(22) 



which leads to an exponential decay on a time 
scale Tp = m/£. Identifying the friction £ with 
1//Z and inserting the mass leads to 



2 R 2 p„ 

TP = 

9 v p v 



(23) 



Now we consider a very short time scale tf, the 
time fluid momentum diffuses a distance 2R, 
i.e. (2R) 2 = 2vt f leading to 



tf 



2R 2 
v 



(24) 



which helps defining the particle Reynolds num- 
ber as 



Re^ 

TS 



Rv, 



(25) 



Finally, we have to discuss another important 
short length scale due to a short range potential 
among the colloidal particles. This scale usually 
determines the maximum time step in Molec- 
ular Dynamics. Guided by the analogy to a 
harmonic oscillator with frequency u> = yk/m, 
we replace the spring constant k with the sec- 
ond derivative of the inter-particle potential 
d 2 V(R)/dR 2 and use the period of this oscil- 
lation to define the interaction time scale 



TV 



2tt 



2tti 



' ml 2 



A 



H 



(26) 



where we approximate the derivative of the po- 
tential by means of the Hamaker-constant Ah 
as a typical size of the potential and a typi- 
cal distance I such as the distance between the 
surface of the particle and the primary poten- 
tial minimum due to the combined effect of van 
der Waals attraction and screened Coulomb re- 
pulsion. Comparison of Ty and Tp can an- 
swer the question, whether the oscillations of 
two particles around the primary or secondary 
minimum are visible or whether the creeping 
or over-damped case is realized where friction 
is dominating over inertia. Analyzing a har- 
monic oscillator with damping constant £ one 
finds that creeping being established at 



Tp < 



Ty 

47r' 



(27) 



In using this relation a lubrication force de- 
scribed in Eq. (@J has to be taken into account. 
This force is proportional to the difference of 
normal velocities of two approaching particles 
and in this sense it can be seen as an additional 
contribution to the friction coefficient. It be- 
comes huge at short inter-particle distances d 
and it will turn out later that even without this 
addition all particles considered here are well 
inside the creeping regime due to the large fric- 
tion in water. This is the justification that so 
far many people used Brownian Dynamics (BD) 
for this system instead of Molecular Dynamics 
In our situation, including ther- 
mal fluctuations and full hydrodynamics con- 
sistently is easier to do in Molecular Dynamics. 
Moreover, with our parameters the MD is at 
least competitive or even faster than previous 
BD calculations. 
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VII. SIMILARITY CONSIDERATIONS 
AND DETERMINATION OF 
SIMULATION PARAMETERS 

A. Introduction 

The determination of parameters for a meso- 
scopic model to quantitatively compare with 
experiment is a non-trivial task. Typical values 
of the parameters in an experiment are listed 
in Tab.Q] For these values of the parameters 
in the experiment all the time scales defined in 
the previous section are calculated and listed in 
Tab.|n] This tells us that the Peclet number is 
Pe = td I ts = 0.74, and we have a competition 
between convection due to gravity and Brown- 
ian motion. The particle Reynolds number is 
very small, i.e. Re = t f /t s = 4.0 • 10~ 7 . The 
ratio of ry to Tp is larger than Air, hence oscilla- 
tions of particles in their short range potentials 
are over-damped, already without considering 
lubrication forces. We get Tp -C tc> since the 
particles are well relaxed before they hit each 
other due to Brownian motion. Tp <C Td, hence 
the transport of momentum through the fluid is 
much faster than if transported directly by the 
particle. These are the dynamical characteris- 
tics which have to be preserved by any param- 
eter changes, in particular, the Peclet-number 
has to be kept exactly the same. Of course, the 
static properties such as the ratio of kinetic en- 
ergy ~ k B T to the potential energies, ~ mgR 
and ~ Ah have to be kept the same too. 
However, using identical parameters as shown 
in Tab.[I]in an MD simulation would require of 
the order of 10rs/rp w 5 • 10 7 iterations to see 
sufficient progress in the sedimentation process. 
This is an unacceptably high numerical effort, 
which must be reduced without significantly 
changing the physics of this process. First, we 
now show how to choose the parameters for a 
simulation using the coupling method I. After 
that, we describe what has to be changed using 
the coupling II. 



B. Determination of the parameters for 
coupling I 

We start by choosing reasonable parame- 
ters for the hydrodynamic part of the code, 
i.e. Stochastic Rotation Dynamics (SRD), since 
this is time-consuming and the most storage- 
intensive part of our simulation. For the mo- 
ment we keep the particle radius constant at 
R = 0.4 /im. Let a be the lattice constant of 
the SRD grid. By choosing a = R/2 a spher- 
ical particle covers about 34 boxes which is a 



sufficient resolution of the particle. We get 
a = 0.2 /im. 

We use an average number of M — 2.5 fluid 
particles per box, which leads to 6M —15 real 
numbers (3 velocity and 3 position coordinates 
in 3 dimensions) to be stored for every box. A 
larger M would reduce Brownian motion and 
increase CPU-time and storage requirements. 
Using a smaller number leads to a very long 
effective mean free path of the fluid particles 
(sometimes there is only one particle per box 
and no collision takes place), which results in a 
large viscosity and a bad resolution of the flow 
field around the colloidal particles. 
Next, we choose the ratio of the mean free path 
A = Tv/fceT '/to/ to the lattice constant, to/ 
is the mass of the fluid particle and T the ef- 
fective temperature of the fluid particles which 
can differ by several orders of magnitude from 
the real temperature of the experiment as will 
be explained later. In Ref. 30] it was discov- 
ered that a ratio A/a smaller than 0.5 leads to 
anomalies in the model, which can be corrected 
by a random shift of the lattice prior to every 
rotation. Here, we set A = 0.6 a = 0.12/im to 
have sufficient resolution of the flow and ran- 
dom shifts are not needed. 
The rotation angle a is taken to be 90° be- 
cause this gives the most simple rotation ma- 
trix. The exact expression for the shear viscos- 
ity for a = 90° is given by [34| 



185t 



1 



1 



-M 



M 



k B T St M + 2 



4m / M 



- r 

(28) 

Inserting M = 2.5 and expressing temperature 
by means of A it follows for our choice of pa- 
rameters that 



v = 0.3052—. 

ot 



(29) 



In order to reproduce the same diffusion coef- 
ficient as seen in experiments, St has to be de- 
termined by means of the Einstein relation, 



D = k B T p = 



k B T 
&-Kvp w R 



(30) 



Setting p w — Mrrif/a 3 , using v from Eq. (|2*9^l 
and expressing k B T /rut by means of A one finds 
St = 0.025a 3 / (DR). Inserting the diffusion co- 
efficient expected in reality from the Einstein 
relation, D = 5.49 • 10~ 13 m 2 /s we arrive at a 
time step St = 0.91 ms for the SRD algorithm. 
This time step is of course too large to resolve 
the motion of colloidal particles due to inter- 
particle forces and friction. Hence, a two-step 
method is needed: The trajectory for the col- 
loidal particles is integrated by another, smaller 
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Particle radius R 


0.4 fim 


Temperature T 


300 K 


mass density of particle p m 


3.9 ■ 10 3 kg/m 3 


mass density of water p w 


1.0 • 10 3 kg/m 3 


Boltzmann constant fce 


1.38 -10 -23 J/K 


kin. viscosity of water v 


10- 6 m 2 /s 


gravity g 


9.81 m/s 2 


Hamaker constant Ah of AI2O3 in H2O 


4.76 ■ 10" 20 J 


distance to primary minimum I 


0.008 /an 



TABLE I: Parameters for the simulation 





Td 


TG 


T V 


TF 


TP 


0.791s 


0.582 s 


49.4 ms 


7.45 [is 


0.320 (j,s 


0.139 ,us 



TABLE II: time scales which arise in a system characterized by the parameters listed in Tab.Q] 



time step StM- This also means, that the ex- 
tensive SRD-procedure is only applied every 
St/StMth iteration of the MD-algorithm, thus 
reducing the required computer power substan- 
tially. 

The way parameters are derived implicitly 
means that we keep rs and td as in reality. 
This corresponds to r s /0.91ms = 869 SRD- 
iterations until a colloidal particle has fallen 
down by one diameter 2R which is affordable. 
The kinematic viscosity in the simulation is 
much smaller than in nature (f mo dei = 1-34 • 
10- n m 2 /s). 

Next, one has to check what happens to the 
particle relaxation time rp. The requirement 
is, that it should be much larger than the one 
given in Tab.[fl (in order to increase numerical 
efficiency) and on the other hand it should still 
be smaller than tq to ensure that particles can 
relax between consecutive collisions caused by 
thermal motion. Following Eq. (|23|) . we obtain 
Tp = 7.69 ms. This is an acceptable value: it 
is much larger than the 0.139 fis seen in reality 
and still smaller than tq — 49.4 ms. Therefore 
it needs 7 SRD-steps to relax a particle which 
means that the process still can be resolved. 
Considering momentum transport in the fluid 
versus direct transport: During time rrj>, mo- 
mentum in the fluid is transported a distance 
x 2 = 2vt g = 33.10a 2 , i.e. x = 5.75a = 2.88R. 
Hence, momentum transport in the fluid is 
only slightly faster than by diffusive transport, 
which is still acceptable, even though in the real 
system it is much faster. This is reflected in a 
Reynolds number which is larger by a factor of 
lO^/^modci = 0.746- 10 5 in the simulation, i.e. 
Re = 2.9 • 10~ 2 . This again reflects the fact 
that the SRD-model is efficient only if Peclet- 
and Reynolds number are in the range between 



0.05 and 20. 

Now, the gravity constant g of the model has to 
be determined requiring that the Stokes veloc- 
ity is the same as given in Tab.Q] Since thermal 
convection of the fluid is not important for our 
simulation, we can neglect gravity on the fluid 
particles. Therefore, there is no buoyancy force 
in the simulation. We can correct for that by as- 
suming a smaller gravity constant modified by 
the density ratio of colloid material and fluid. 
We find 

k'modcl / Pw \ 

.9modol = ffrcal 1 I 

freal \ Pm J 

= 9.78 -lO" 5 -^ (31) 
S^ 

As mentioned above, not only the viscosity, but 
also the temperature in our simulation may be 
different from the one in nature. To see that 
we calculate the ratio A = Pm/k^T . In nature 
we have A = 0.942 • 10 24 s 2 /m 5 . In the model 
we get Amodci = 3.9 Mrrif/(a 3 k^T) where we 
express ksT/mf by means of the mean free 
path and the time step A 2 /(St) 2 . One finds 
that A mo doi is scaled by a factor of 7.44 • 10 4 . 
The static features have to be reproduced by 
the model, and therefore we have to keep the 
ratio of kinetic and potential energy k^T/An 
constant. This means the ratio p m /Vp ot and 
especially p m /Aji has also to be scaled by this 
factor. We use A H = 4.76-10 - 20 J/(7.44-10 4 ) in 
the model, corresponding to new Ah = 8.61 • 
10~ 25 J. From Eq. f2^jl we get a scaled ry of 
2.03 ms corresponding to Ty/rp = 0.264 (which 
is smaller than 47r, see Ea. l27|) . The unsealed 
value is 53.6. The creeping case is restored by 
the lubrication force, which we have included 
in the MD simulation, and which grows for 
smaller gaps between the particles. The lubri- 
cation force determines the small iteration time 



10 



step StM for the MD simulation. We chose 
5tM = 2/is, which is about 200 times larger 
than it would be if all the original parameters 
would have been kept and mhi(Ty, rp), being 
much smaller, would determine the time step. 
Comparing to the SRD time step we see that ev- 
ery 455 small steps one SRD step is performed. 
We need 869 SRD steps and 4.3 • 10 6 MD steps 
to see a colloidal particle sinking down by one 
diameter. The time scales in the simulation are 
summarized again in Tab.lml 

C. Determination of the parameters for 
coupling II 

To simulate the same system with coupling 
method II, we use the same particle radius 
R = 0.4 /j to. The lattice constant has now to 
be chosen differently because the colloidal parti- 
cles are coupled to the SRD-simulation as mass 
points. They have influence on the fluid which 
is in the same cell, and therefore the size of the 
cell can be understood as the volume within 
which the SRD-simulation "feels" the colloidal 
particles and we choose the lattice constant in 
a way that the volume of the cell is equal to the 
volume of a colloidal particle: a — 6.25 • 10~ 7 m. 
A smaller lattice constant in this context would 
model smaller colloidal particles in the SRD- 
part of the simulation. The velocity field would 
be resolved better, but since coupling method 
II does not allow a resolution smaller than the 
colloidal particles, one can not expect to gain 
any information from the fluid simulation on 
smaller length scales than the colloidal particle 
size. Any attempt to increase the resolution of 
the SRD simulation would only cause a larger 
computational effort. 

Since we do not modify the Peclet number, we 
have to choose approximately the same number 
of fluid particles per colloidal particle. Since the 
box size has increased with respect to the cou- 
pling method I, we have to assume more parti- 
cles per box now. We choose M — 60 (which 
would correspond to two particles per box in 
the coupling method I, but since the boxes are 
much larger now, we can slightly reduce the ra- 
tio of fluid particles per colloidal particle). 
We choose A/a = 0.5 and use random grid shifts 
here to avoid that fluid particles interact too of- 
ten with the same partners which causes arte- 
facts in their correlation. The rotation angle 
a is again 90° to achieve very simple matrices. 
Following the same procedure as for coupling 
method I fEa. l28fl we find a time step for the 
SRD of 

St = 2.05 ms. (32) 



According to Eq.|2Uthe viscosity in the simula- 
tion results to ^ mo dci = 2.29 • 10~ n m 2 /s. The 
gravity constant therefore has to be rescaled by 
a factor of 58813, and the temperature and the 
potentials have to be scaled by 43733 43] . 

The resulting characteristic times are shown 
in Tab. II VI t$ and td are again kept as in re- 
ality. Now we need r s /0.00106s = 385 SRD- 
iterations until a colloidal particle has fallen 
down by one diameter, which is much faster 
than by using coupling method I. Tp is still 
smaller than tq, but these two times are now of 
the same order of magnitude. This reflects the 
fact that with coupling method II the local flow 
around the particles cannot be resolved. Relax- 
ation of the particles between their collisions 
due to Brownian motion is in this case less im- 
portant because lubrication effects can not be 
seen here, tf is of comparable size, too, which 
means that the diffusion of the momentum is 
now on the same time scale as the particle mo- 
tion. This is understandable since the length 
scales of the colloidal particles and of the fluid 
boxes have to be the same. 

Momentum is transported 1.2 times faster in 
the fluid as by the particles themselves. Short 
range hydrodynamic interactions which cannot 
be resolved are in this sense comparable to par- 
ticle - particle collisions whereas for long range 
interactions the slightly faster transport of mo- 
mentum can reproduce coarse grained hydro- 
dynamic effects. Again, to model these effects 
comparable to reality, the Reynolds number 
has to be much smaller than unity. We find 
Re = 1.77- 10- 2 . 

If we include lubrication forces in the MD 
simulation in order to reproduce at least to 
some extend short range hydrodynamics, we 
have to choose the same MD time step as for 
coupling method I but we need approximately 
fifty per cent less CPU time for the hydrody- 
namics. Even though it seems to be a too 
simplified approach, we can reproduce a vol- 
ume fraction dependent sedimentation velocity 
as will be described in the results section. 



VIII. SIMULATION SETUP 

A. Boundary conditions 

Most simulations have been performed using 
periodic boundary conditions in all three di- 
rections. Then the total momentum may not 
change in any simulation step if no external 
forces (like gravity) are applied. If gravity on 
the colloidal particles is applied in a system 
with periodic boundary conditions, this would 
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1~S 


TD 




TV 


TF 


Tp 


0.791 s 


0.582 s 


49.4 ms 


2.03 ms 


22.9 ms 


7.69 ms 



TABLE III: Time scales in the simulation using coupling method I. 



TS 


Td 






TF 


Tp 


0.791s 


0.582 s 


49.4 ms 


1.56 ms 


14.0 ms 


18.2 ms 



TABLE IV: Time scales in the simulation using coupling method II. 



accelerate the whole system, since the total 
force on the center of mass is not vanishing. In 
a real system there is friction at the walls and 
even more important, there is an equilibrium 
between hydrostatic pressure acting on the sur- 
face of a given volume and the gravity acting as 
a body force. Since we simulate a volume in the 
center of the suspension we either have to ap- 
ply the pressure on the walls or, which is easier, 
make sure that in sum the forces on the center 
of mass of the whole simulated system vanishes. 
Therefore, we follow the center of mass, i.e. on 
particles with higher density their gravity mi- 
nus buoyancy has to be applied, so that they 
move downward whereas the same force in op- 
posite direction has to be applied to the fluid, 
which makes it move upward like in a sedimen- 
tation vessel with a closed bottom. 

For the following discussion we define that 
the direction in which eventually gravity is ap- 
plied is called — z-direction, if a shear force is 
applied acts in the x-direction. Using closed 
boundaries wall effects may be introduced, e.g. 
crystallization starts earlier than in the bulk. 
This effect could be observed especially when 
gravity was switched off and only closed bound- 
ary conditions were applied. This is a finite 
size effect, which is not that strong, if periodic 
boundaries are applied. But in the case of grav- 
ity being applied, the whole system accelerates. 
To face this problem, three possibilities were 
tested: 

1) fix the boundaries only in z-direction, 

2) fix the boundaries in x and y-direction and 
apply no-slip for the fluid, 

3) choose periodic boundaries in all directions 
and compensate the gravitation on the colloidal 
particles with a force in the opposite direction 
applied on the fluid. 

Possibilities 1) and 2) simulate a system close 
to a wall, in case 1) it is the bottom of a ves- 
sel whereas in case 2) the experiment would be 
done in a capillary. Possibility 3) turned out to 
be the most realistic simulation although it can 
start to drift, if the compensating force is not 
adjusted accurately. Slowly accumulated drifts 
of the center of mass can be removed every hun- 



dreds of SRD-time steps if necessary. 



B. Temperature and thermostat 

We have measured the temperature of the 
colloidal particles for different setups If damp- 
ing constants are chosen appropriately, the re- 
sulting temperature fits very well the tempera- 
ture, which we have adjusted for the fluid by the 
initial conditions. If we additionally switch on 
a thermostat which we describe in the follow- 
ing the measured temperature exactly agrees 
with temperature adjusted by the thermostat. 
When gravity is applied to the system, parti- 
cles are accelerated and if in addition periodic 
boundaries are used, a thermostat is absolutely 
needed to remove the extra energy, introduced 
by the periodic boundary in z direction in com- 
bination with gravity. 

Therefore we use a modified version of 
the thermostat described in [2(j [Chap. 7.4.1, 
"Stochastic methods" p. 227 f.]. The thermo- 
stat, originally suggested by Heves|39|. chooses 
a random scaling factor £ for the velocities from 
an interval [1 — 7,1 + 7]. The scaling of the ve- 
locity is then accepted or rejected according to 
a Monte Carlo scheme. However, the detailed 
balance is not fulfilled for the choice of ( de- 
scribed in 26]. In our implementation of the 
thermostat, we randomly choose an e in the in- 
terval [0, 7] and apply for £ one of the values 
1 + e or , each of them with the probabil- 
ity of \. With one of these values the velocity 
is scaled by the Monte Carlo acceptance rate. 
Also the temperature in our case is defined 
slightly different from[2^]: the mean velocity 
u within one SRD-cell defines the velocity field 
of the fluid and gives the hydrodynamic inter- 
action between the colloidal particles. There- 
fore it may not be modified by the thermostat. 
We only scale the velocity component relative 
to the mean velocity: v" 6 " 1 = £(vj — u) + u. 
The Monte Carlo acceptance rate in our case is 
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given by 

C (3(M-1)) exp (_( M _ !)(£2 _ !) T / T *) 

M (W) 

r= 5p ^E(v i -u)^, 

i— 1 

which is the local temperature in the SRD-cell 
and T* denotes the temperature to which the 
thermostat will drive the system. M is the 
number of particles in the cell. Note that one 
has to divide the total thermal energy in the 
SRD-cell by M— 1 instead of M to calculate the 
local temperature. This reflects the fact that 
the mean velocity u in the cell already contains 
three degrees of freedom which the particles in 
the SRD-cell have. The choice of 7 and the fre- 
quency with which the thermostat is called to 
work determine the relaxation rate, with which 
the system adapts T*. The version described 
in [2^| shows deviations of the achieved temper- 
ature for small numbers of particles per cell, 
whereas our implementation exactly reproduces 
T*. The thermostat can even be extended to 
particles of different mass i.e. colloid and fluid 
particles where the mass is used as weight fac- 
tor for all velocities of the simulation. 

C. Outlook: Shear 

There are several possibilities to shear the 
system. If one only has MD particles, one can 
use moving walls either with a spring constant 
and a friction coefficient or with direct hard re- 
flections, where a moving wall is assumed and 
the reflection is calculated in the moving frame 
fixed to the wall. 

These approaches of course neglect all effects 
(like pseudo wall slip), which appear close to a 
wall in a shear experiment with a suspension. 
There, shear stress has to be applied to the fluid 
which then drags the suspended particles. One 
way to implement this, is to add a small ve- 
locity offset to all fluid particles which are re- 
flected. Since this approach works well and the 
colloidal particles are dragged by the fluid, we 
apply shear in this way to our system. 

IX. TESTS OF THE SIMULATION 
CODE 

A. Conservation of energy, velocity 
distributions 

We have checked that the total energy is con- 
served in the molecular dynamics simulation if 
all damping constants are switched off. Oth- 
erwise, or if the total energy even increases in 



spite of damping constants, the MD time step 
has been chosen too large. In the SRD simula- 
tion energy is conserved as well and if we use 
coupling method II also for the total system 
energy is conserved within numerical accuracy. 
With coupling method I (where a thermostat is 
already included in the coupling method) or if 
we switch on an additional thermostat energy 
will not exactly be conserved but the system 
will reach a stable, i.e. equilibrated state. In 
that sense, total energy (including thermal en- 
ergy) will converge to a constant value. 
In SRD-Simulations without any embedded 
particles, the total energy contains only the ki- 
netic energy of the fluid particles. It is fully 
determined by the initialization of the particle 
velocities. We can choose three uniformly dis- 
tributed random numbers to initialize the three 
velocity components for the fluid particles. In 
thermal equilibrium the distribution should be 
a Gaussian, which in fact can be observed in 
our simulations after some tens of SRD time 
steps. If colloidal particles are included into 
the system, they should reach a thermal equi- 
librium, at least as long as no external forces 
are applied. Damping terms would reduce fluc- 
tuations, so, to check, if the colloidal particles 
reach the same temperature as the fluid parti- 
cles, damping constants have to be set to zero. 
Both distributions are shown in Fig. [21 They 
are both Gaussian with the correct tempera- 
ture, even though for initialization uniformly 
distributed random numbers (square well) had 
been used. The tests are performed with both 
coupling methods. We have carried out simula- 
tions with particle radii of 0.4 fim. and 0.25 /im, 
where the Peclet number (for the simulations 
where gravity is applied) is 0.11 and of course, 
it takes much longer to observe sedimentation. 



B. Viscosity 

The diffusion coefficient of suspended col- 
loidal particles can be used to check if the de- 
sired viscosity could really be achieved in the 
simulation. Using Eq. I|18|) we can, once we have 
measured D, calculate the kinematic viscosity 
v and compare it to the value we have used to 
determine the simulation parameters like the 
SRD time step. We achieve a deviation of less 
than 20% in a diluted system compared to the 
theoretical value for an infinitely diluted sys- 
tem. Note that D is a fixed number only in the 
limit of an infinitely diluted system and only if 
the interaction potentials between the colloidal 
particles are exclusively repulsive. 
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FIG. 2: Velocity distribution of fluid (a) and colloid 
(b) particles in a SRD simulation after thermaliza- 
tion. Particle density is 3900 kg/m 3 , the time step 
is 2.0 /xs, model temperature T — 10.57mK, fluid 
particle mass m/ = 1.0667 ■ 10 -18 , particle diame- 
ter d = 0.5 jj,m. The theoretical Gaussian curve is 
plotted as well as the measured velocity distribu- 
tions. 



We are using two different methods, either 
the Green-Kubo-method or direct evaluation of 
the mean square displacement. The first is even 
very accurate, if only few particles are used, but 
consumes much computer time and memory be- 
cause all particle velocities have to be stored 
for all time steps used in the calculation. That 
means, for higher volume fractions, it is more 
efficient just to sum up all the mean square dis- 
placements within a given period of time. To 
calculate D using the Green-Kubo method one 
uses the following relation: 



9x{j) 



lim — — (34) 
/^oo IM To t 

/ M T ot 

^ ^ v x . n ((i + j)St)v Xtn (i5t) 

i—1 n— 1 



D, 



St f^x(0) + £sx0')j 



(35) 



where Mx t is the total number of particles in 
the system, / is the number of time steps used 



to calculate the contribution g(j). v x , n (i5t) de- 
notes the x component of the velocity of particle 
n in the i-th time step. The sum in the expres- 
sion for D x is in principle an infinite one, but 
since the contributions g(j) decay with j~ 3 ^ 2 , 
one can truncate this sum after some tens of 
terms. D y and D z can be calculated accord- 
ingly. In Fig. 0a we show the diffusion coef- 
ficient in each direction. In numerical calcu- 
lations it is impossible to evaluate an infinite 
sum. In Eq. Ij34(l / is limited at least by the to- 
tal number of time steps within the simulation 
and in Eq. i|35fl the sum therefore is not infinite 
either. Since the contributions g x (j) become 
more and more inaccurate for larger j we trun- 
cate the sum after n terms and find that in our 
simulations for n « 50 the diffusion coefficient 
does not change anymore in a systematic way 
if n is increased further. In Fig. 0b the last 
term of the sum is shown. For larger values, 
they fluctuate due to the finite sum in Eq. (|34|l 
which leads to the inaccuracy in the right part 
of Fig.[3]a. These fluctuations become smaller 
for longer simulation runs, but do not change 
the value of the diffusion coefficient taken as an 
average from the center part of Fig. 0a 44] . 
For the mean square displacement in one direc- 
tion during a time interval At we calculate 



D x = 



1 



A/Tot 



2AtM Tot 



(x l (t + At) - x t (t)f 



(36) 

and D y and D z accordingly. For medium den- 
sities we have compared both methods and 
achieved the same results within error bars. De- 
pending on the number of particles, we use one 
of both methods. 

According to Richardson and Zaki 40J , the 
mean sedimentation velocity of particles sus- 
pended in a liquid depends on the volume frac- 
tion <b as: 



>(i 



(37) 



with a typical exponent I between « 2.5 and 
4 depending on the boundary conditions. For 
periodic boundary conditions, Peclet number of 
Pe = 1 and Reynolds number Re <§C 1 we find 
an exponent of 3.5 (Fig.0J even when we use 
coupling method II, where only long range hy- 
drodynamic interaction can be calculated cor- 
rectly. A similar value is found for Pe = 2 
and Pe = \. Padding and Louis have found 
that the exponent I depends very weakly on the 
Peclet number |23|. We have used here coupling 
method II, but some investigations have also 
been carried out using coupling method I. On 
the first view there is no big difference apparent 
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a) Evaluation of D using the Green-Kubo method: the plot shows the sum of 



IJj=i. 



and the estimated D 



b) the decay of the contributions g x , y ,z(j)- We have measured the diffusion constant of soft 
spheres coupled with coupling II to the SRD at low volume fractions. 
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FIG. 4: Mean sedimentation velocity over porosity 
(1 — 0) according to Eq. (1371) : Measured values and 
fit curve in a log- log-plot. The Peclet number of 
this simulation is 1. Coupling method II has been 
used for this plot. 



between the two coupling methods, at least, as 
long as, like in this test of our simulation code, 
no attractive forces are included. Our first re- 
sults where we have studied the peloid system 
in more detail are presented in the following 
section. 



X. RESULTS 

A. Spacial correlation functions 

For our production runs, we have simulated vol- 
ume fractions of 7, 14, 21, 28 and 35% in a 
cube with an extension of (6/im) 3 . Therefore 
are 231, 462, 693, 924 and 1155 colloidal par- 
ticles respectively and 2.0 • 10 5 fluid particles 
necessary. 

We have evaluated the particle-particle corre- 



lation function. For attractive potentials sev- 
eral sharp peaks can be observed and we assign 
them to distinct local orders of particles. Os- 
cillations can be found in the correlation func- 
tion. They are caused by exclusion of volume. 
In the case of attractive forces they are less pro- 
nounced than if mainly repulsive interaction is 
present, see Fig. |5] and Fig.|H| 
In Fig. [5] the particles cluster due to their at- 
tractive potentials and form stable configura- 
tions. The diameter of the particles is 5T0~ 7 m. 
There is a sharp peak in the spacial correla- 
tion function of the particle centers at exactly 
that distance 2R, where two particles touch 
each other in the very left part of the plot (A). 
Then, for larger particle separation, the correla- 
tion function starts to grow and drops suddenly 
after a peak at 1 /im (D), which is twice the di- 
ameter (4-R). This is the contribution of two 
particles touching the same third particle. The 
distance between them depends on the angle, 
which they form with the particle in the mid- 
dle, but, it is at last twice the diameter, when 
they are in a straight line, which explains the 
sudden drop of the correlation function. If sev- 
eral particles stick together, the straight line is 
stabilized. This explains the peak at the end of 
this section of the correlation function. 
Two more peaks can clearly be assigned to con- 
figurations: One of them is from two particles 
touching two other particles, which themselves 
touch each other (C). There again the case of 
all particles being in the same plane can be 
stabilized by other particles surrounding them. 
The particles under consideration are then sep- 
arated by a distance of 2i?v / 3- But of course, 
bending this configuration is still a degree of 
freedom which brings the two particles slightly 
closer to each other. Thus their contribution to 



15 




LJ ■ ■ ' ' ■ ■ ■ ■ 1 

2345 6789 10 11 

center - center distance / particle radius 

FIG. 5: Correlation function of AI2O3 for $q = 50 mV and k — 3 ■ 10 8 /m. The potential is attractive, thus 
peaks (labeled by letters) can be identified and assigned to special local configurations (see text). 
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FIG. 6: Correlation function of AI2O3 for ^0 = 50 mV and k — 7.3 ■ 10 7 /m. repulsive potentials. One can 
see oscillations caused by the excluded volume. 



the correlation function is shifted downward. 
The fourth peak at reflects two particles, 

both touching three particles, which themselves 
are touching each other and define a plane (B). 
There is no freedom anymore for the two parti- 
cles touching all the three of them at the same 
time. One can place one of them at one side of 



the plane and the other one at the other side. 
When the potentials are mainly repulsive and 
the minimum caused by the van der Waals at- 
traction is only a fraction of fceT, the spatial 
correlation function looks completely different, 
as depicted in Fig.0 The peaks described in 
the previous paragraphs have disappeared here. 
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FIG. 7: Plot of the correlation function of Fig.E] 
together with the potential used in this simulation. 
One can see that the maximum of the correlation 
function occurs for the distance, at which the very 
shallow secondary minimum of the potential is lo- 
cated. 



The primary peak has moved to a slightly larger 
distance, since the repulsive potential hinders 
the particles from touching each other. 
In Fig.JJJwe compare the correlation function of 
Fig.[S]with the potential used for that simula- 
tion. The maximum of the correlation function 
coincides with the minimum of the potential, 
but, as the minimum is not very sharp, the par- 
ticles are not restricted to fixed geometries and 
are in a steady process of rearrangement which 
results in broader peaks. This process could 
also be studied by evaluating the velocity corre- 
lation function for the colloidal particles which 
is related to the viscosity of the sample. The 
correlation of particles which arc several diam- 
eters apart is still remarkable, as it is transmit- 
ted by the particles in between. The oscillations 
of the correlation function can be understood as 
a formation of layers where the probability of 
finding a particle in a certain layer is higher 
than in between. 



B. Shear 

We have carried out simulations with shear 
and gravity. For the particles the boundaries in 
z direction were closed, gravity was applied in 
negative ^-direction only to the colloidal par- 
ticles. For the fluid particles the boundary in 
z-direction was closed as well and additionally 
a velocity offset was added to apply a shear in 
x-direction. Boundaries for fluid and for parti- 
cles were periodic in x- and y-direction. Veloc- 
ity distribution functions have been evaluated. 
For the cases we investigated, after a transient 
they are all Gaussian (Fig.|SJ). 
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FIG. 8: Velocity distribution of colloidal particles 
for each direction. Semi-log-plot where deviations 
from a Gaussian would be visible by deviations 
from a parabolic profile. 



C. Phase diagram 

We have explored the phase diagram for 
AI2O3 with respect to screening length and ef- 
fective surface potential. We could identify the 
regions of suspended single particles and of fioc- 
culation (Fig.[5|). The transition between these 
two regions depends on both parameters, Debye 
screening length and effective surface potential. 
It is known that the pH- value determines the ef- 
fective surface potential 'J'o, and that salt con- 
centration and pH-value determine the Debye 
screening length k . Exact relations between 
salt concentration and pH- value on one side and 
k and on the other side are not known a pri- 
ori for the parameter ranges of our suspensions. 
There are approximations for very diluted sys- 
tems and low salt concentrations. It is known 
that for AI2O3 the surface potential becomes 
zero for pH « 8.7^l]. However, a phase transi- 
tion between clustering in the upper left part of 
Fig.lHland a suspended regime in the lower right 
part can be found in the simulations in analogy 
to the experiment. The spatial correlation func- 
tion can be evaluated for all the simulated cases 
and it can be used as a tool to identify the two 
regions of the phase diagram. 
Figs. llUtTHfl show selected examples of correla- 
tion functions for different parameter sets. The 
first and second graph refer to a volume frac- 
tion $ = 14% which also has been used for the 
phase diagram of Fig.EU In Fig.^Hthe correla- 
tion function has been plotted for every other 
image of the left column in the phase diagram 
in Fig.0 One can see that for suspended par- 
ticles only the first peak can be found in the 
correlation function. The secondary minimum 
in the potential causes the particles to glue for 
short times before they continue with their dif- 
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30 50 V rff (mV) 

FIG. 9: Snapshots from the phase diagram of AI2O3: For the DLVO-Potentials with different effective 
surface charge and different screening length one can either observe cluster formation or single particles in 
suspension. The simulation was done at room temperature for 1 second of real time and a particle diameter 
of 0.5 fj,m. Gravity has not been applied here. The pictures are corresponding to the values written on 
the axis. For this figure we have chosen the simulation runs for 14% volume fraction with 462 colloidal 
particles. 



fusion process. With increasing k the secondary 
minimum approaches the particle surface, and 
therefore the main peak is shifted to smaller 
distances. At the same time it becomes deeper 
so that clusters are formed and more peaks oc- 
cur. The peak at a distance of 2R^/2 1=3 3i? 
disappears again, when the attraction becomes 
stronger since this is a meta stable configura- 
tion of particles forming an octahedron. Fig. llll 
corresponds to the first row of images of Fig. [3] 
In this case the depth of the secondary mini- 
mum is adjusted by changing the effective sur- 
face potential. Again the transition between 
clustering regime and suspension can be ob- 
served. The potentials used here are among the 



ones plotted above in Fig.^45]. In Fig. 1 121 and 
IT31 the dependence of the correlation function 
on the volume fraction can be seen. In both 
cases long range correlations become more pro- 
nounced with increasing volume fraction. This 
is shown for the suspended regime fFig. 11211 and 
for the clustering regime fFig.H3[l. where the 
transition between the two cases presented here 
is achieved by a variation of n by only 10%. 



D. Diffusion 

We measured the diffusion coefficient of col- 
loidal particles with attractive potentials. In 
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FIG. 10: Correlation function and its dependence 
on the inverse Debye screening length k. = 
20 mV and "!> = 0.14 have been kept constant. For 
shorter Debye screening lengths the attractive force 
becomes stronger and leads to clustering, which is 
reflected in the appearance of peaks. The single 
curves have been shifted with respect to each other. 



FIG. 12: Correlation function and its dependence 
on the volume fraction $. Effective surface poten- 
tial * = 20 mV and k = 1.4 ■ 10 8 m" 1 have been 
kept constant. For center-center distances between 
six and eight particle radii broad peaks start to ap- 
pear for larger volume fractions. 
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FIG. 11: Correlation function and its depen- 
dence on the effective surface potential ^o- it = 
2 ■ 10 s m _1 and $ = 0.14 have been kept con- 
stant. The higher the effective surface potential, 
the stronger the attraction force and clustering can 
be seen in the growing peaks. 
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FIG. 13: Correlation function and its dependence 
on the volume fraction $. Effective surface poten- 
tial * = 20 mV and k = 1.6 ■ 10 8 m" 1 have been 
kept constant. Due to a small change in n with 
respect to Fig. 1121 one can cross the phase border 
between suspended particles and clustering regime. 
Also here long range correlations become more pro- 
nounced for high volume fractions. 



Fig-El we show the diffusion coefficient for 
AI2O3 with an effective surface potential of 
= 50 mV and an inverse Debye screening 
length of k = 2 • 10 8 m -1 for room temperature. 
One can see that the mobility of the particles 
decays since a cluster formation process takes 
place and the particles in the cluster are rela- 
tively fixed. The remaining mobility consists of 
two parts: Particles can still, with a non vanish- 
ing probability, leave the cluster by thermal ac- 
tivation and the cluster itself can take part in a 
diffusion process, it can vibrate or be deformed 
- all of these are processes which are taking 
place on much longer time scales than the single 
particle diffusion. By studying the dependency 
of the diffusion coefficient on the potentials and 



on the volume fraction, one might be able to 
find an answer to the question, which of these 
processes is important for the dynamics of the 
system in which part of the phase diagram of 
Fig.® 



XI. CONCLUSION 

We have shown that by combining a Stochastic 
Rotation Dynamics and a Molecular Dynamics 
simulation it is possible to study dense colloidal 
suspensions. We have explained how to de- 
termine effective parameters for the simulation 
(box size a, simulation time step St, number 
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FIG. 14: Mean square displacement (projection on 
each of the axis). For this simulation $o = 50 mV 
and k = 3 • 10 8 /m have been used. First parti- 
cles move by diffusion, are attracted and then form 
clusters of larger size and lower mobility, which can 
be observed in a decay of the diffusion length for a 
given period of time. We have simulated 1155 par- 
ticles in a cube with 6 pm extension, which results 
in a volume fraction of 35%. 



of fluid particles per box M . . .). It is possible 
to relate the simulation to very distinct exper- 
imental conditions since all parameters (den- 
sity, temperature, potentials. . . ) which enter 
into the description are scaled in a well defined 
manner. We have presented first results which 
demonstrate the power of the model. We have 
demonstrated that the Richardson-Zaki law is 
reproduced already with the simple and fast 



coupling method II and we have studied the 
dependence of the pair correlation function on 
the shape of the interaction potentials. We have 
shown how one can distinguish if for given De- 
bye screening length k, effective surface poten- 
tial \&o an d Hamaker constant an if the system 
is in the clustering or suspended regime. 
We are planning to carry out detailed investi- 
gations of the properties described in the two 
preceding sections (diffusion coefficient, corre- 
lation functions, sedimentation velocity) as well 
as cluster size and shape. Then these quantities 
can be analyzed under shear, their dependence 
on the shear rate, and the shear viscosity of the 
suspension, containing the fluid and the parti- 
cles, which both contribute to a complex shear 
viscosity. 
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